home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 16 / CU Amiga Magazine's Super CD-ROM 16 (1997-10-16)(EMAP Images)(GB)[!][issue 1997-11].iso / CUCD / Graphics / Ghostscript / source / gsmisc.c < prev    next >
C/C++ Source or Header  |  1997-06-10  |  24KB  |  816 lines
  1. /* Copyright (C) 1989, 1996, 1997 Aladdin Enterprises.  All rights reserved.
  2.   
  3.   This file is part of Aladdin Ghostscript.
  4.   
  5.   Aladdin Ghostscript is distributed with NO WARRANTY OF ANY KIND.  No author
  6.   or distributor accepts any responsibility for the consequences of using it,
  7.   or for whether it serves any particular purpose or works at all, unless he
  8.   or she says so in writing.  Refer to the Aladdin Ghostscript Free Public
  9.   License (the "License") for full details.
  10.   
  11.   Every copy of Aladdin Ghostscript must include a copy of the License,
  12.   normally in a plain ASCII text file named PUBLIC.  The License grants you
  13.   the right to copy, modify and redistribute Aladdin Ghostscript, but only
  14.   under certain conditions described in the License.  Among other things, the
  15.   License requires that the copyright notice and this notice be preserved on
  16.   all copies.
  17. */
  18.  
  19. /* gsmisc.c */
  20. /* Miscellaneous utilities for Ghostscript library */
  21. #include "malloc_.h"
  22. #include "math_.h"
  23. #include "memory_.h"
  24. #include "gx.h"
  25. #include "gpcheck.h"        /* for gs_return_check_interrupt */
  26. #include "gserrors.h"
  27. #include "gconfigv.h"        /* for USE_ASM */
  28. #include "gxfarith.h"
  29. #include "gxfixed.h"
  30.  
  31. /* Define private replacements for stdin, stdout, and stderr. */
  32. FILE *gs_stdin, *gs_stdout, *gs_stderr;
  33. /* Ghostscript writes debugging output to gs_debug_out. */
  34. /* We define gs_debug and gs_debug_out even if DEBUG isn't defined, */
  35. /* so that we can compile individual modules with DEBUG set. */
  36. char gs_debug[128];
  37. FILE *gs_debug_out;
  38.  
  39. /* Define eprintf_program_name and lprintf_file_and_line as procedures */
  40. /* so one can set breakpoints on them. */
  41. void
  42. eprintf_program_name(FILE *f, const char *program_name)
  43. {    fprintf(f, "%s: ", program_name);
  44. }
  45. void
  46. lprintf_file_and_line(FILE *f, const char *file, int line)
  47. {    fprintf(f, "%s(%d): ", file, line);
  48. }
  49.  
  50. /* Log an error return.  We always include this, in case other */
  51. /* modules were compiled with DEBUG set. */
  52. #undef gs_log_error        /* in case DEBUG isn't set */
  53. int
  54. gs_log_error(int err, const char _ds *file, int line)
  55. {    if ( gs_log_errors )
  56.       { if ( file == NULL )
  57.           dprintf1("Returning error %d.\n", err);
  58.         else
  59.           dprintf3("%s(%d): Returning error %d.\n",
  60.                (const char *)file, line, err);
  61.       }
  62.     return err;
  63. }
  64.  
  65. /* Check for interrupts before a return. */
  66. int
  67. gs_return_check_interrupt(int code)
  68. {    if ( code < 0 )
  69.       return code;
  70.     { int icode = gp_check_interrupts();
  71.       return (icode == 0 ? 0 :
  72.           gs_note_error((icode > 0 ? gs_error_interrupt : icode)));
  73.     }
  74. }
  75.  
  76. /* ------ Substitutes for missing C library functions ------ */
  77.  
  78. #ifdef memory__need_memmove        /* see memory_.h */
  79. /* Copy bytes like memcpy, guaranteed to handle overlap correctly. */
  80. /* ANSI C defines the returned value as being the src argument, */
  81. /* but with the const restriction removed! */
  82. void *
  83. gs_memmove(void *dest, const void *src, size_t len)
  84. {    if ( !len )
  85.         return (void *)src;
  86. #define bdest ((byte *)dest)
  87. #define bsrc ((const byte *)src)
  88.     /* We use len-1 for comparisons because adding len */
  89.     /* might produce an offset overflow on segmented systems. */
  90.     if ( ptr_le(bdest, bsrc) )
  91.       {    register byte *end = bdest + (len - 1);
  92.         if ( ptr_le(bsrc, end) )
  93.           {    /* Source overlaps destination from above. */
  94.             register const byte *from = bsrc;
  95.             register byte *to = bdest;
  96.             for ( ; ; )
  97.               {    *to = *from;
  98.                 if ( to >= end )    /* faster than = */
  99.                   return (void *)src;
  100.                 to++; from++;
  101.               }
  102.           }
  103.       }
  104.     else
  105.       {    register const byte *from = bsrc + (len - 1);
  106.         if ( ptr_le(bdest, from) )
  107.           {    /* Source overlaps destination from below. */
  108.             register const byte *end = bsrc;
  109.             register byte *to = bdest + (len - 1);
  110.             for ( ; ; )
  111.               {    *to = *from;
  112.                 if ( from <= end )    /* faster than = */
  113.                   return (void *)src;
  114.                 to--; from--;
  115.               }
  116.           }
  117.       }
  118. #undef bdest
  119. #undef bsrc
  120.     /* No overlap, it's safe to use memcpy. */
  121.     memcpy(dest, src, len);
  122.     return (void *)src;
  123. }
  124. #endif
  125.  
  126. #ifdef memory__need_memchr        /* see memory_.h */
  127. /* ch should obviously be char rather than int, */
  128. /* but the ANSI standard declaration uses int. */
  129. const char *
  130. gs_memchr(const char *ptr, int ch, size_t len)
  131. {    if ( len > 0 )
  132.     {    register const char *p = ptr;
  133.         register uint count = len;
  134.         do { if ( *p == (char)ch ) return p; p++; } while ( --count );
  135.     }
  136.     return 0;
  137. }
  138. #endif
  139.  
  140. #ifdef memory__need_memset        /* see memory_.h */
  141. /* ch should obviously be char rather than int, */
  142. /* but the ANSI standard declaration uses int. */
  143. void *
  144. gs_memset(void *dest, register int ch, size_t len)
  145. {    if ( ch == 0 )
  146.         bzero(dest, len);
  147.     else if ( len > 0 )
  148.     {    register char *p = dest;
  149.         register uint count = len;
  150.         do { *p++ = (char)ch; } while ( --count );
  151.     }
  152.     return dest;
  153. }
  154. #endif
  155.  
  156. #ifdef malloc__need_realloc        /* see malloc_.h */
  157. /* Some systems have non-working implementations of realloc. */
  158. void *
  159. gs_realloc(void *old_ptr, size_t old_size, size_t new_size)
  160. {    void *new_ptr;
  161.     if ( new_size )
  162.       { new_ptr = malloc(new_size);
  163.         if ( new_ptr == NULL )
  164.           return NULL;
  165.       }
  166.     else
  167.       new_ptr = NULL;
  168.     /* We have to pass in the old size, since we have no way to */
  169.     /* determine it otherwise. */
  170.     if ( old_ptr != NULL )
  171.       { if ( new_ptr != NULL )
  172.           memcpy(new_ptr, old_ptr, min(old_size, new_size));
  173.         free(old_ptr);
  174.       }
  175.     return new_ptr;
  176. }
  177. #endif
  178.  
  179. /* ------ Debugging support ------ */
  180.  
  181. /* Dump a region of memory. */
  182. void
  183. debug_dump_bytes(const byte *from, const byte *to, const char *msg)
  184. {    const byte *p = from;
  185.     if ( from < to && msg )
  186.         dprintf1("%s:\n", msg);
  187.     while ( p != to )
  188.        {    const byte *q = min(p + 16, to);
  189.         dprintf1("0x%lx:", (ulong)p);
  190.         while ( p != q ) dprintf1(" %02x", *p++);
  191.         dputc('\n');
  192.        }
  193. }
  194.  
  195. /* Dump a bitmap. */
  196. void
  197. debug_dump_bitmap(const byte *bits, uint raster, uint height, const char *msg)
  198. {    uint y;
  199.     const byte *data = bits;
  200.     for ( y = 0; y < height; ++y, data += raster )
  201.         debug_dump_bytes(data, data + raster, (y == 0 ? msg : NULL));
  202. }
  203.  
  204. /* Print a string. */
  205. void
  206. debug_print_string(const byte *chrs, uint len)
  207. {    uint i;
  208.     for ( i = 0; i < len; i++ )
  209.       dputc(chrs[i]);
  210.     fflush(dstderr);
  211. }
  212.  
  213. /* ------ Arithmetic ------ */
  214.  
  215. /* Compute M modulo N.  Requires N > 0; guarantees 0 <= imod(M,N) < N, */
  216. /* regardless of the whims of the % operator for negative operands. */
  217. int
  218. imod(int m, int n)
  219. {    if ( n <= 0 )
  220.       return 0;        /* sanity check */
  221.     if ( m >= 0 )
  222.       return m % n;
  223.     { int r = -m % n;
  224.       return (r == 0 ? 0 : n - r);
  225.     }
  226. }
  227.  
  228. /* Compute the GCD of two integers. */
  229. int
  230. igcd(int x, int y)
  231. {    int c = x, d = y;
  232.     if ( c < 0 ) c = -c;
  233.     if ( d < 0 ) d = -d;
  234.     while ( c != 0 && d != 0 )
  235.       if ( c > d ) c %= d;
  236.       else d %= c;
  237.     return d + c;        /* at most one is non-zero */
  238. }
  239.  
  240. #if defined(set_fmul2fixed_vars) && !USE_ASM
  241.  
  242. /*
  243.  * Floating multiply with fixed result, for avoiding floating point in
  244.  * common coordinate transformations.  Assumes IEEE representation,
  245.  * 16-bit short, 32-bit long.  Optimized for the case where the first
  246.  * operand has no more than 16 mantissa bits, e.g., where it is a user space
  247.  * coordinate (which are often integers).
  248.  *
  249.  * The assembly language version of this code is actually faster than
  250.  * the FPU, if the code is compiled with FPU_TYPE=0 (which requires taking
  251.  * a trap on every FPU operation).  If there is no FPU, the assembly
  252.  * language version of this code is over 10 times as fast as the emulated FPU.
  253.  */
  254. /* Some of the following code has been tweaked for the Borland 16-bit */
  255. /* compiler.  The tweaks do not change the algorithms. */
  256. #if arch_ints_are_short && !defined(FOR80386)
  257. #  define SHORT_ARITH
  258. #endif
  259. int
  260. set_fmul2fixed_(fixed *pr, long /*float*/ a, long /*float*/ b)
  261. {
  262. #ifdef SHORT_ARITH
  263. #  define long_rsh8_ushort(x)\
  264.     (((ushort)(x) >> 8) | ((ushort)((ulong)(x) >> 16) << 8))
  265. #  define utemp ushort
  266. #else
  267. #  define long_rsh8_ushort(x) ((ushort)((x) >> 8))
  268. #  define utemp ulong
  269. #endif
  270.     /* utemp may be either ushort or ulong.  This is OK because */
  271.     /* we only use ma and mb in multiplications involving */
  272.     /* a long or ulong operand. */
  273.     utemp ma = long_rsh8_ushort(a) | 0x8000;
  274.     utemp mb = long_rsh8_ushort(b) | 0x8000;
  275.     int e = 260 + _fixed_shift - ((
  276.         (((uint)((ulong)a >> 16)) & 0x7f80) +
  277.         (((uint)((ulong)b >> 16)) & 0x7f80)
  278.       ) >> 7);
  279.     ulong p1 = ma * (b & 0xff);
  280.     ulong p = (ulong)ma * mb;
  281. #define p_bits (size_of(p) * 8)
  282.  
  283.     if ( (byte)a )        /* >16 mantissa bits */
  284.     {    ulong p2 = (a & 0xff) * mb;
  285.         p += ((((uint)(byte)a * (uint)(byte)b) >> 8) + p1 + p2) >> 8;
  286.     }
  287.     else
  288.         p += p1 >> 8;
  289.     if ( (uint)e < p_bits )        /* e = -1 is possible */
  290.         p >>= e;
  291.     else if ( e >= p_bits )        /* also detects a=0 or b=0 */
  292.       {    *pr = fixed_0;
  293.         return 0;
  294.       }
  295.     else if ( e >= -(p_bits - 1) || p >= 1L << (p_bits - 1 + e) )
  296.         return_error(gs_error_limitcheck);
  297.     else
  298.         p <<= -e;
  299.     *pr = ((a ^ b) < 0 ? -p : p);
  300.     return 0;
  301. }
  302. int
  303. set_dfmul2fixed_(fixed *pr, ulong /*double lo*/ xalo, long /*float*/ b, long /*double hi*/ xahi)
  304. {
  305. #ifdef SHORT_ARITH
  306. #  define long_lsh3(x) ((((x) << 1) << 1) << 1)
  307. #  define long_rsh(x,ng16) ((uint)((x) >> 16) >> (ng16 - 16))
  308. #else
  309. #  define long_lsh3(x) ((x) << 3)
  310. #  define long_rsh(x,ng16) ((x) >> ng16)
  311. #endif
  312.     return set_fmul2fixed_(pr,
  313.                    (xahi & 0xc0000000) +
  314.                     (long_lsh3(xahi) & 0x3ffffff8) +
  315.                     long_rsh(xalo, 29),
  316.                    b);
  317. }
  318.  
  319. #endif
  320.  
  321. #if USE_FPU_FIXED
  322.  
  323. /*
  324.  * Convert from floating point to fixed point with scaling.
  325.  * These are efficient algorithms for FPU-less machines.
  326.  */
  327. #define mbits_float 23
  328. #define mbits_double 20
  329. int
  330. set_float2fixed_(fixed *pr, long /*float*/ vf, int frac_bits)
  331. {    fixed mantissa;
  332.     int shift;
  333.  
  334.     if ( !(vf & 0x7f800000) )
  335.       { *pr = fixed_0;
  336.         return 0;
  337.       }
  338.     mantissa = (fixed)((vf & 0x7fffff) | 0x800000);
  339.     shift = ((vf >> 23) & 255) - (127 + 23) + frac_bits;
  340.     if ( shift >= 0 )
  341.     {    if ( shift >= sizeof(fixed) * 8 - 24 )
  342.           return_error(gs_error_limitcheck);
  343.         if ( vf < 0 )
  344.           mantissa = -mantissa;
  345.         *pr = (fixed)(mantissa << shift);
  346.     }
  347.     else
  348.       *pr = (shift < -24 ? fixed_0 :
  349.          vf < 0 ? -(fixed)(mantissa >> -shift) :  /* truncate */
  350.          (fixed)(mantissa >> -shift));
  351.     return 0;
  352. }
  353. int
  354. set_double2fixed_(fixed *pr, ulong /*double lo*/ lo,
  355.   long /*double hi*/ hi, int frac_bits)
  356. {    fixed mantissa;
  357.     int shift;
  358.  
  359.     if ( !(hi & 0x7ff00000) )
  360.       { *pr = fixed_0;
  361.         return 0;
  362.       }
  363.     /* We only use 31 bits of mantissa even if sizeof(long) > 4. */
  364.     mantissa = (fixed)(((hi & 0xfffff) << 10) | (lo >> 22) | 0x40000000);
  365.     shift = ((hi >> 20) & 2047) - (1023 + 30) + frac_bits;
  366.     if ( shift > 0 )
  367.       return_error(gs_error_limitcheck);
  368.     *pr = (shift < -30 ? fixed_0 :
  369.            hi < 0 ? -(fixed)(mantissa >> -shift) :  /* truncate */
  370.            (fixed)(mantissa >> -shift));
  371.     return 0;
  372. }
  373. /*
  374.  * Given a fixed value x with fbits bits of fraction, set v to the mantissa
  375.  * (left-justified in 32 bits) and f to the exponent word of the
  376.  * corresponding floating-point value with mbits bits of mantissa in the
  377.  * first word.  (The exponent part of f is biased by -1, because we add the
  378.  * top 1-bit of the mantissa to it.)
  379.  */
  380. static const byte f2f_shifts[] =
  381.  { 4, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0 };
  382. #define f2f_declare(v, f)\
  383.     ulong v;\
  384.     long f
  385. #define f2f(x, v, f, mbits, fbits)\
  386.     if ( x < 0 )\
  387.       f = 0xc0000000 + (29 << mbits) - ((long)fbits << mbits), v = -x;\
  388.     else\
  389.       f = 0x40000000 + (29 << mbits) - ((long)fbits << mbits), v = x;\
  390.     if ( v < 0x8000 )\
  391.       v <<= 15, f -= 15 << mbits;\
  392.     if ( v < 0x800000 )\
  393.       v <<= 8, f -= 8 << mbits;\
  394.     if ( v < 0x8000000 )\
  395.       v <<= 4, f -= 4 << mbits;\
  396.     { int shift = f2f_shifts[v >> 28];\
  397.       v <<= shift, f -= shift << mbits;\
  398.     }
  399. long
  400. fixed2float_(fixed x, int frac_bits)
  401. {    f2f_declare(v, f);
  402.  
  403.     if ( x == 0 )
  404.       return 0;
  405.     f2f(x, v, f, mbits_float, frac_bits);
  406.     return f + (((v >> 7) + 1) >> 1);
  407. }
  408. void
  409. set_fixed2double_(double *pd, fixed x, int frac_bits)
  410. {    f2f_declare(v, f);
  411.  
  412.     if ( x == 0 )
  413.       { ((long *)pd)[1 - arch_is_big_endian] = 0;
  414.         ((ulong *)pd)[arch_is_big_endian] = 0;
  415.       }
  416.     else
  417.       { f2f(x, v, f, mbits_double, frac_bits);
  418.         ((long *)pd)[1 - arch_is_big_endian] = f + (v >> 11);
  419.         ((ulong *)pd)[arch_is_big_endian] = v << 21;
  420.       }
  421. }
  422.  
  423. /*
  424.  * Compute A * B / C when 0 <= B < C and A * B exceeds (or might exceed)
  425.  * the capacity of a long.
  426.  */
  427. #ifdef DEBUG
  428. struct { long mnanb, mnab, manb, mab, mnc, mdq, mde, mds, mqh, mql; } fmq_stat;
  429. #  define mincr(x) ++fmq_stat.x
  430. #else
  431. #  define mincr(x) DO_NOTHING
  432. #endif
  433. fixed
  434. fixed_mult_quo(fixed signed_A, fixed B, fixed C)
  435. {    /* First compute A * B in double-fixed precision. */
  436.     ulong A = (signed_A < 0 ? -signed_A : signed_A);
  437.     long msw;
  438.     ulong lsw;
  439.     ulong p1;
  440.  
  441. #define num_bits (sizeof(fixed) * 8)
  442. #define half_bits (num_bits / 2)
  443. #define half_mask ((1L << half_bits) - 1)
  444.     if ( B <= half_mask )
  445.       { if ( A <= half_mask )
  446.           { fixed Q = (ulong)(A * B) / (ulong)C;
  447.  
  448.             mincr(mnanb);
  449.         return (signed_A < 0 ? -Q : Q);
  450.           }
  451.         /*
  452.          * We might still have C <= half_mask, which we can
  453.          * handle with a simpler algorithm.
  454.          */
  455.         lsw = (A & half_mask) * B;
  456.         p1 = (A >> half_bits) * B;
  457.         if ( C <= half_mask )
  458.           { ulong q0 = (p1 += lsw >> half_bits) / C;
  459.             ulong rem = ((p1 - C * q0) << half_bits) + (lsw & half_mask);
  460.         ulong Q = (q0 << half_bits) + rem / C;
  461.  
  462.         mincr(mnc);
  463.         return (signed_A < 0 ? -Q : Q);
  464.           }
  465.         msw = p1 >> half_bits;
  466.         mincr(manb);
  467.       }
  468.     else if ( A <= half_mask )
  469.       { p1 = A * (B >> half_bits);
  470.         msw = p1 >> half_bits;
  471.         lsw = A * (B & half_mask);
  472.         mincr(mnab);
  473.       }
  474.     else
  475.       { /* We have to compute all 4 products.  :-( */
  476.         ulong lo_A = A & half_mask;
  477.         ulong hi_A = A >> half_bits;
  478.         ulong lo_B = B & half_mask;
  479.         ulong hi_B = B >> half_bits;
  480.         ulong p1x = hi_A * lo_B;
  481.  
  482.         msw = hi_A * hi_B;
  483.         lsw = lo_A * lo_B;
  484.         p1 = lo_A * hi_B;
  485.         if ( p1 > max_ulong - p1x )
  486.           msw += 1L << half_bits;
  487.         p1 += p1x;
  488.         msw += p1 >> half_bits;
  489.         mincr(mab);
  490.       }
  491.     /* Finish up by adding the low half of p1 to the high half of lsw. */
  492. #if max_fixed < max_long
  493.     p1 &= half_mask;
  494. #endif
  495.     p1 <<= half_bits;
  496.     if ( p1 > max_ulong - lsw )
  497.       msw++;
  498.     lsw += p1;            
  499.     /*
  500.      * Now divide the double-length product by C.  Note that we know msw
  501.      * < C (otherwise the quotient would overflow).  Start by shifting
  502.      * (msw,lsw) and C left until C >= 1 << (num_bits - 1).
  503.      */
  504.     { ulong denom = C;
  505.       int shift = 0;
  506.  
  507. #define bits_4th (num_bits / 4)
  508.       if ( denom < 1L << (num_bits - bits_4th) )
  509.         { mincr(mdq);
  510.           denom <<= bits_4th, shift += bits_4th;
  511.         }
  512. #undef bits_4th
  513. #define bits_8th (num_bits / 8)
  514.       if ( denom < 1L << (num_bits - bits_8th) )
  515.         { mincr(mde);
  516.           denom <<= bits_8th, shift += bits_8th;
  517.         }
  518. #undef bits_8th
  519.       while ( !(denom & (1L << (num_bits - 1))) )
  520.         { mincr(mds);
  521.           denom <<= 1, ++shift;
  522.         }
  523.       msw = (msw << shift) + (lsw >> (num_bits - shift));
  524.       lsw <<= shift;
  525. #if max_fixed < max_long
  526.       lsw &= (1L << (sizeof(fixed) * 8)) - 1;
  527. #endif
  528.       /* Compute a trial upper-half quotient. */
  529.       { ulong hi_D = denom >> half_bits;
  530.         ulong lo_D = denom & half_mask;
  531.         ulong hi_Q = (ulong)msw / hi_D;
  532.         /* hi_Q might be too high by 1 or 2, but it isn't too low. */
  533.         ulong p0 = hi_Q * hi_D;
  534.         ulong p1 = hi_Q * lo_D;
  535.         ulong hi_P;
  536.  
  537.         while ( (hi_P = p0 + (p1 >> half_bits)) > msw ||
  538.             (hi_P == msw && ((p1 & half_mask) << half_bits) > lsw)
  539.           )
  540.           { /* hi_Q was too high by 1. */
  541.         --hi_Q;
  542.         p0 -= hi_D;
  543.         p1 -= lo_D;
  544.         mincr(mqh);
  545.           }
  546.         p1 = (p1 & half_mask) << half_bits;
  547.         if ( p1 > lsw )
  548.           msw--;
  549.         lsw -= p1;
  550.         msw -= hi_P;
  551.         /* Now repeat to get the lower-half quotient. */
  552.         msw = (msw << half_bits) + (lsw >> half_bits);
  553. #if max_fixed < max_long
  554.         lsw &= half_mask;
  555. #endif
  556.         lsw <<= half_bits;
  557.         { ulong lo_Q = (ulong)msw / hi_D;
  558.           long Q;
  559.  
  560.           p1 = lo_Q * lo_D;
  561.           p0 = lo_Q * hi_D;
  562.           while ( (hi_P = p0 + (p1 >> half_bits)) > msw ||
  563.               (hi_P == msw && ((p1 & half_mask) << half_bits) > lsw)
  564.             )
  565.         { /* lo_Q was too high by 1. */
  566.           --lo_Q;
  567.           p0 -= hi_D;
  568.           p1 -= lo_D;
  569.           mincr(mql);
  570.         }
  571.           Q = (hi_Q << half_bits) + lo_Q;
  572.           return (signed_A < 0 ? -Q : Q);
  573.         }
  574.       }
  575.     }
  576. #undef half_bits
  577. #undef half_mask
  578. }
  579.  
  580. #endif
  581.  
  582. /* Trace calls on sqrt when debugging. */
  583. #undef sqrt
  584. extern double sqrt(P1(double));
  585. double
  586. gs_sqrt(double x, const char *file, int line)
  587. {    if ( gs_debug_c('~') )
  588.       {    fprintf(stdout, "[~]sqrt(%g) at %s:%d\n",
  589.             x, (const char *)file, line);
  590.         fflush(stdout);
  591.       }
  592.     return sqrt(x);
  593. }
  594.  
  595. /*
  596.  * Define sine and cosine functions that take angles in degrees rather than
  597.  * radians, and that are implemented efficiently on machines with slow
  598.  * (or no) floating point.
  599.  */
  600. #if USE_FPU < 0            /****** maybe should be <= 0 ? ******/
  601.  
  602. #define sin0 0.00000000000000000
  603. #define sin1 0.01745240643728351
  604. #define sin2 0.03489949670250097
  605. #define sin3 0.05233595624294383
  606. #define sin4 0.06975647374412530
  607. #define sin5 0.08715574274765817
  608. #define sin6 0.10452846326765346
  609. #define sin7 0.12186934340514748
  610. #define sin8 0.13917310096006544
  611. #define sin9 0.15643446504023087
  612. #define sin10 0.17364817766693033
  613. #define sin11 0.19080899537654480
  614. #define sin12 0.20791169081775931
  615. #define sin13 0.22495105434386498
  616. #define sin14 0.24192189559966773
  617. #define sin15 0.25881904510252074
  618. #define sin16 0.27563735581699916
  619. #define sin17 0.29237170472273671
  620. #define sin18 0.30901699437494740
  621. #define sin19 0.32556815445715670
  622. #define sin20 0.34202014332566871
  623. #define sin21 0.35836794954530027
  624. #define sin22 0.37460659341591201
  625. #define sin23 0.39073112848927377
  626. #define sin24 0.40673664307580015
  627. #define sin25 0.42261826174069944
  628. #define sin26 0.43837114678907740
  629. #define sin27 0.45399049973954675
  630. #define sin28 0.46947156278589081
  631. #define sin29 0.48480962024633706
  632. #define sin30 0.50000000000000000
  633. #define sin31 0.51503807491005416
  634. #define sin32 0.52991926423320490
  635. #define sin33 0.54463903501502708
  636. #define sin34 0.55919290347074679
  637. #define sin35 0.57357643635104605
  638. #define sin36 0.58778525229247314
  639. #define sin37 0.60181502315204827
  640. #define sin38 0.61566147532565829
  641. #define sin39 0.62932039104983739
  642. #define sin40 0.64278760968653925
  643. #define sin41 0.65605902899050728
  644. #define sin42 0.66913060635885824
  645. #define sin43 0.68199836006249848
  646. #define sin44 0.69465837045899725
  647. #define sin45 0.70710678118654746
  648. #define sin46 0.71933980033865108
  649. #define sin47 0.73135370161917046
  650. #define sin48 0.74314482547739413
  651. #define sin49 0.75470958022277201
  652. #define sin50 0.76604444311897801
  653. #define sin51 0.77714596145697090
  654. #define sin52 0.78801075360672190
  655. #define sin53 0.79863551004729283
  656. #define sin54 0.80901699437494745
  657. #define sin55 0.81915204428899180
  658. #define sin56 0.82903757255504174
  659. #define sin57 0.83867056794542394
  660. #define sin58 0.84804809615642596
  661. #define sin59 0.85716730070211222
  662. #define sin60 0.86602540378443860
  663. #define sin61 0.87461970713939574
  664. #define sin62 0.88294759285892688
  665. #define sin63 0.89100652418836779
  666. #define sin64 0.89879404629916704
  667. #define sin65 0.90630778703664994
  668. #define sin66 0.91354545764260087
  669. #define sin67 0.92050485345244037
  670. #define sin68 0.92718385456678731
  671. #define sin69 0.93358042649720174
  672. #define sin70 0.93969262078590832
  673. #define sin71 0.94551857559931674
  674. #define sin72 0.95105651629515353
  675. #define sin73 0.95630475596303544
  676. #define sin74 0.96126169593831889
  677. #define sin75 0.96592582628906831
  678. #define sin76 0.97029572627599647
  679. #define sin77 0.97437006478523525
  680. #define sin78 0.97814760073380558
  681. #define sin79 0.98162718344766398
  682. #define sin80 0.98480775301220802
  683. #define sin81 0.98768834059513777
  684. #define sin82 0.99026806874157036
  685. #define sin83 0.99254615164132198
  686. #define sin84 0.99452189536827329
  687. #define sin85 0.99619469809174555
  688. #define sin86 0.99756405025982420
  689. #define sin87 0.99862953475457383
  690. #define sin88 0.99939082701909576
  691. #define sin89 0.99984769515639127
  692. #define sin90 1.00000000000000000
  693.  
  694. private const double sin_table[361] = {
  695.   sin0,
  696.   sin1, sin2, sin3, sin4, sin5, sin6, sin7, sin8, sin9, sin10,
  697.   sin11, sin12, sin13, sin14, sin15, sin16, sin17, sin18, sin19, sin20,
  698.   sin21, sin22, sin23, sin24, sin25, sin26, sin27, sin28, sin29, sin30,
  699.   sin31, sin32, sin33, sin34, sin35, sin36, sin37, sin38, sin39, sin40,
  700.   sin41, sin42, sin43, sin44, sin45, sin46, sin47, sin48, sin49, sin50,
  701.   sin51, sin52, sin53, sin54, sin55, sin56, sin57, sin58, sin59, sin60,
  702.   sin61, sin62, sin63, sin64, sin65, sin66, sin67, sin68, sin69, sin70,
  703.   sin71, sin72, sin73, sin74, sin75, sin76, sin77, sin78, sin79, sin80,
  704.   sin81, sin82, sin83, sin84, sin85, sin86, sin87, sin88, sin89, sin90,
  705.   sin89, sin88, sin87, sin86, sin85, sin84, sin83, sin82, sin81, sin80,
  706.   sin79, sin78, sin77, sin76, sin75, sin74, sin73, sin72, sin71, sin70,
  707.   sin69, sin68, sin67, sin66, sin65, sin64, sin63, sin62, sin61, sin60,
  708.   sin59, sin58, sin57, sin56, sin55, sin54, sin53, sin52, sin51, sin50,
  709.   sin49, sin48, sin47, sin46, sin45, sin44, sin43, sin42, sin41, sin40,
  710.   sin39, sin38, sin37, sin36, sin35, sin34, sin33, sin32, sin31, sin30,
  711.   sin29, sin28, sin27, sin26, sin25, sin24, sin23, sin22, sin21, sin20,
  712.   sin19, sin18, sin17, sin16, sin15, sin14, sin13, sin12, sin11, sin10,
  713.   sin9, sin8, sin7, sin6, sin5, sin4, sin3, sin2, sin1, sin0,
  714.   -sin1, -sin2, -sin3, -sin4, -sin5, -sin6, -sin7, -sin8, -sin9, -sin10,
  715. -sin11, -sin12, -sin13, -sin14, -sin15, -sin16, -sin17, -sin18, -sin19, -sin20,
  716. -sin21, -sin22, -sin23, -sin24, -sin25, -sin26, -sin27, -sin28, -sin29, -sin30,
  717. -sin31, -sin32, -sin33, -sin34, -sin35, -sin36, -sin37, -sin38, -sin39, -sin40,
  718. -sin41, -sin42, -sin43, -sin44, -sin45, -sin46, -sin47, -sin48, -sin49, -sin50,
  719. -sin51, -sin52, -sin53, -sin54, -sin55, -sin56, -sin57, -sin58, -sin59, -sin60,
  720. -sin61, -sin62, -sin63, -sin64, -sin65, -sin66, -sin67, -sin68, -sin69, -sin70,
  721. -sin71, -sin72, -sin73, -sin74, -sin75, -sin76, -sin77, -sin78, -sin79, -sin80,
  722. -sin81, -sin82, -sin83, -sin84, -sin85, -sin86, -sin87, -sin88, -sin89, -sin90,
  723. -sin89, -sin88, -sin87, -sin86, -sin85, -sin84, -sin83, -sin82, -sin81, -sin80,
  724. -sin79, -sin78, -sin77, -sin76, -sin75, -sin74, -sin73, -sin72, -sin71, -sin70,
  725. -sin69, -sin68, -sin67, -sin66, -sin65, -sin64, -sin63, -sin62, -sin61, -sin60,
  726. -sin59, -sin58, -sin57, -sin56, -sin55, -sin54, -sin53, -sin52, -sin51, -sin50,
  727. -sin49, -sin48, -sin47, -sin46, -sin45, -sin44, -sin43, -sin42, -sin41, -sin40,
  728. -sin39, -sin38, -sin37, -sin36, -sin35, -sin34, -sin33, -sin32, -sin31, -sin30,
  729. -sin29, -sin28, -sin27, -sin26, -sin25, -sin24, -sin23, -sin22, -sin21, -sin20,
  730. -sin19, -sin18, -sin17, -sin16, -sin15, -sin14, -sin13, -sin12, -sin11, -sin10,
  731.   -sin9, -sin8, -sin7, -sin6, -sin5, -sin4, -sin3, -sin2, -sin1, -sin0
  732. };
  733.  
  734. double
  735. gs_sin_degrees(double ang)
  736. {    int ipart;
  737.     if ( is_fneg(ang) )
  738.       ang = 180 - ang;
  739.     ipart = (int)ang;
  740.     if ( ipart >= 360 )
  741.       { int arem = ipart % 360;
  742.         ang -= (ipart - arem);
  743.         ipart = arem;
  744.       }
  745.     return
  746.       (ang == ipart ? sin_table[ipart] :
  747.        sin_table[ipart] + (sin_table[ipart+1] - sin_table[ipart]) *
  748.          (ang - ipart));
  749. }
  750.  
  751. double
  752. gs_cos_degrees(double ang)
  753. {    int ipart;
  754.     if ( is_fneg(ang) )
  755.       ang = 90 - ang;
  756.     else
  757.       ang += 90;
  758.     ipart = (int)ang;
  759.     if ( ipart >= 360 )
  760.       { int arem = ipart % 360;
  761.         ang -= (ipart - arem);
  762.         ipart = arem;
  763.       }
  764.     return
  765.       (ang == ipart ? sin_table[ipart] :
  766.        sin_table[ipart] + (sin_table[ipart+1] - sin_table[ipart]) *
  767.          (ang - ipart));
  768. }
  769.  
  770. void
  771. gs_sincos_degrees(double ang, gs_sincos_t *psincos)
  772. {    psincos->sin = gs_sin_degrees(ang);
  773.     psincos->cos = gs_cos_degrees(ang);
  774. }
  775.  
  776. #else                /* we have floating point */
  777.  
  778. static const int isincos[5] = { 0, 1, 0, -1, 0 };
  779.  
  780. double
  781. gs_sin_degrees(double ang)
  782. {    double quot = ang / 90;
  783.     if ( floor(quot) == quot )
  784.       { int quads = (int)fmod(quot, 4) & 3;    /* & 3 because might be < 0 */
  785.         return isincos[quads];
  786.       }
  787.     return sin(ang * (M_PI / 180));
  788. }
  789.  
  790. double
  791. gs_cos_degrees(double ang)
  792. {    double quot = ang / 90;
  793.     if ( floor(quot) == quot )
  794.       { int quads = (int)fmod(quot, 4) & 3;    /* & 3 because might be < 0 */
  795.         return isincos[quads + 1];
  796.       }
  797.     return cos(ang * (M_PI / 180));
  798. }
  799.  
  800. void
  801. gs_sincos_degrees(double ang, gs_sincos_t *psincos)
  802. {    double quot = ang / 90;
  803.     if ( floor(quot) == quot )
  804.       { int quads = (int)fmod(quot, 4) & 3;    /* & 3 because might be < 0 */
  805.         psincos->sin = isincos[quads];
  806.         psincos->cos = isincos[quads + 1];
  807.       }
  808.     else
  809.       { double arad = ang * (M_PI / 180);
  810.         psincos->sin = sin(arad);
  811.         psincos->cos = cos(arad);
  812.       }
  813. }
  814.  
  815. #endif                /* USE_FPU */
  816.